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demiological models with ordinary derivatives replaced by fractional deriva¬ 
tives in an an-hoc manner. These models may be mathematically interesting 
but their relevance is uncertain. Here we develop an SIR model for an epidemic, 
including vital dynamics, from an underlying stochastic process. We show how 
fractional differential operators arise naturally in these models whenever the 
recovery time from the disease is power law distributed. This can provide a 
model for a chronic disease process where individuals who are infected for a 
long time are unlikely to recover. The fractional order recovery model is shown 
to be consistent with the Kermack-McKendrick age-structured SIR model and 
it reduces to the Hethcote-Tudor integral equation SIR model. The deriva¬ 
tion from a stochastic process is extended to discrete time, providing a stable 
numerical method for solving the model equations. 

We have carried out simulations of the fractional order recovery model 
showing convergence to equilibrium states. The number of infecteds in the 
endemic equilibrium state increases as the fractional order of the derivative 
tends to zero. 
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1 Introduction 

The classic SIR model was introduced by Kermack and McKendrick in 1927 
[2'5j . In the simplest form of this model the population is separated into three 
compartments representing susceptible (S), infected (I) and removed (R) pop¬ 
ulations. The time evolution of the disease is then described by three coupled 
ordinary differential equations. Kermack and McKendrick also considered more 
general age-structured SIR models in 1932 and 1933 1261127]. These more gen¬ 
eral models were formulated as a set of coupled integro-differential equations. 
Numerous other variants of the SIR model, and related models, have been 
widely studied over the past several decades [32]. The classical SIR models 
have also been formulated as the expectation of stochastic processes jj. 

It is easy to find numerous studies in the literature of ad-hoc fractional or¬ 
der SIR models where the time derivatives are simply replaced with nonlocal 
Caputo fractional derivatives EMlStlllllliaEESlElESlIia ■ These models, 
and their analysis, are mathematically interesting, however they have no par¬ 
ticular physical motivation, other than the incorporation of a memory effect. 
Spatially dependent SIR models with space fractional derivatives have also 
been considered [19] to account for nonlocal spread of disease. 

In this work we show that fractional order time derivatives can be included 
in the governing equations for an SIR model derived from a physical underly¬ 
ing stochastic process. We consider a variation of an SIR model where recovery 
from the disease is dependent on the time since infection. The model is derived 
from a directed continuous time random walk through the SIR compartments, 
with the time in the infectious compartment drawn from a waiting time prob¬ 
ability density. We show that, in the case of a power law tailed waiting time 
density, the governing equations become a set of fractional order differential 
equations. The expected recovery time diverges in a power law waiting time 
density and this leads to chronic infection in the fractional order SIR model. 
As the fractional order derivative operates on the recovery we refer to this as 
the frSIR model. Other fractional epidemic compartment models, such as a 
fractional SEIR model, could also be derived following the approach outlined 
here. There have been several studies of semi-Markovian epidemic models in 
the recent literature ESGH that are related to the approach presented here, 
but they are not formulated as coupled integro-differential equations. 

Starting with a discrete time stochastic process formulation of the frSIR 
model we derive a numerical method for solving the governing fractional or¬ 
der differential equations. The numerical method is related to the discrete 
time stochastic process method that was recently introduced to solve the frac¬ 
tional Fokker-Planck equation [3] . We have implemented the numerical scheme 
to investigate the effects of changes in the fractional order exponent on the 
qualitative behaviour of solutions. The numerical solutions converge to the 
calculated equilibrium states of the frSIR model, when the parameters are 
constants. 

In section [2] we derive an SIR model with an arbitrary waiting time before 
transitioning from the infectious compartment to the recovered compartment. 
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This general model is shown to be consistent with the structured formulation 
of Kernrack and McKendrick [26lf27] . We also derive an integral equation rep¬ 
resentation of the model and show that it reduces to the integral equations 
presented by Hethcote and Tudor [2T] when the parameters are constants. In 
section [3] we show that in the case of a power-law waiting time distribution 
in the infectious compartment we obtain fractional order derivatives in the 
model and we present the governing equations for the frSIR model. In section 
[4] a stable numerical scheme for solving the frSIR model is derived from a dis¬ 
crete time formulation of the stochastic process. In Appendix B the discrete 
time formulation is shown to converge, under a continuous time limit, to the 
continuous time formulation. Numerical solutions are investigated in section 
[bj and we conclude with a summary in section [6] 


2 SIR as a Continuous Time Random Walk 

An underlying assumption in the simplest SIR models is that the transition of 
an individual through each of the compartments is independent of the amount 
of time since the individual entered the compartment. This assumption is 
mathematically equivalent to assuming that the time spent in each compart¬ 
ment is exponentially distributed. This is a very restrictive assumption with 
no a priori reason for it to hold. In some diseases with the potential for chronic 
infection, such as human papillomavirus (HPV), there is evidence of power- 
law tails in the distribution of infected times We have incorporated an 
arbitrary time in the infected compartment in our derivation of a generalised 
SIR model below by way of a continuous time random walk (CTRW) [H2113£j . 

In general our derivation may be adapted to any compartment model where 
the transition out of a compartment is dependent on the length of time since 
entering the compartment. In this work we have concentrated on an SIR model 
with births and deaths. 

In the standard manner, we separate the population into three compart¬ 
ments, Susceptible (S), Infectious (I), and Recovered (R) J5]. The population 
is composed of individuals who are born into the S compartment and undergo 
a directed CTRW on the S, I, and R, compartments until they die and are 
removed from consideration. As in the standard model, individuals may only 
move from the S compartment to the I compartment and then to the R com¬ 
partment. The transition to the I compartment occurs when an individual 
becomes infected and the transition to the R compartment occurs when an 
individual recovers from the infection. The derivation of fractional diffusion 
equations mm, fractional reaction diffusion equations mmm, fractional 
Fokker-Planck equations mmm and fractional cliemotaxis diffusion equations 
Eng from CTRWs has been well studied and provides clear physical moti¬ 
vation for each of these systems. Our derivation of the evolution equations for 
the SIR model and fractional SIR model from a stochastic CTRW process with 
reactions is similar to the derivation of the fractional Fokker-Planck equation 
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with reactions [2] and the derivation of the master equations for CTRWs with 
reactions on networks [3]. 

Consider an individual who is infectious. The probability that they will 
infect a susceptible person in the time t to t + St is assumed to be a product 
of the probability that the infectious person will encounter a susceptible and 
the probability that an encounter with a susceptible will result in an infection. 
Without loss of generality we can express the probability that an encounter 
with a susceptible will result in infection by ui(t)St + o(St ), identifying uj(t) as 
the rate of becoming infected per time interval St. Given that there are S(t) 
susceptible people at time t, this implies that the probability of an infected 
individual creating a new infected individual in the time interval t to t + St is 
uj(t)S(t)St + o(St). We represent the flux of individuals entering the infected 
state at time t, by q + (I, t), which can be recursively constructed from the flux 
at earlier times. Explicitly we have 



u}(t)S(t)@(t, t')q + (I , t')dt', 


(1) 


where d>{t : t') is the probability that an infected individual has survived in the 
infected state until time t given that they entered the state at time t!. Let 
i(—t 0) be the number of individuals who became infected at time t! < 0 and 
who are still infected at time 0, hence, 




<P(0,t’) ’ 


t' < 0. 


( 2 ) 


We can then write Eq. (JT|) for t > 0 as, 

q+(I,t)=[ uj{t)S(t)^{t,t/)q + (I,t')dt' + [ v{t)S(t) ^ i(-t ', 0 )d£. 

JO J-oo 1 ) 

(3) 

For an individual to be in the infected compartment at time t they must 
have become infected at time t or at some prior time t! and remained in the 
compartment until t. The number infected at time t can therefore be found 
from the flux, and the survival probability, as follows, 


I(t) = I 0 (t) + 


<P(t, t')q + (I , t')dt'. 


(4) 


Here we have defined the function, 


Io(t) 



g (*,0 

<P(0,t') 


i(—t', 0 )dt'. 


(5) 


We assume that there are two possible ways in which an individual can move 
from the infectious compartment; they can either recover from the disease and 
move to the R compartment, or they can die and be removed from considera¬ 
tion. If these two possibilities are independent we may write, 


, t') = (j)(t — t')9(t , t') 


( 6 ) 
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where (f>(t — t') is the probability of surviving the jump transition to the R 
compartment from time t! to time t, and 0(0 t ') is the probability of surviving 
death from time t' until time t. If we have a time dependent death rate such 
that the probability of death occurring in the interval t to t+St is 7 (t)St + o(8t) 
then the death survival can be written as, 

0{t,t') = e-R'*W da . (7) 

The rate of change of the population in the infectious compartment can be 
found by differentiating Eq. Q to get, 

= q+ ( I , t ) ~ J 'ip(t-t , )e(t,t , )q + (I,t')dt' - 7 (i) J (!>(t-t')0{t,t')q + {I,t')dt' + < ^^~ 

= g + (40 - i/j{t -t')9(t,t')q + {i,t')dt' - 7(040 + e{t,o)^ ■ 

( 8 ) 


Here we have used the fact that, </>(0) = 1, and that the derivative of the 
jump survival function, <p(t), is a waiting time probability density function, 
here denoted i.e., 


d<j>(t) 

dt 




(9) 


Substituting Eq. <§ into Eq. ^ gives, 


dl(t) 

dt 


=u(t)S(t) (/ 4>(t - t')q + {I, t')dt' + Io(t)j 

I ^ t ~ t '') e ( t ' t '') q+ ( I ' t ') dt ' - 7(040 + 0 ( 00 )^ 


4(0 

0(t, 0 ) 


( 10 ) 


In order to obtain a generalised master equation we need to express the right 
hand side of this equation in terms of I(t). We first use the definition of I(t) 
in Eq. Q, to write, 


dI ^ } - = u(t)S(t)I{t)~ J ^(t-t , )9(t,t')q + (I,t/)dt'-'y(t)I(t)+9(t,0)^ 

Further, noting that, 


4(0 \ 
0(t,0)J ■ 
( 11 ) 


0(0 0) = 0(0 t')9(t', 0) V0 <t' <t, 


we can write Eq. © as, 


m 

0(0 0 ) 


4(0 

0 ( 00 ) 



L /xff + ( 40) 
J 0(o,o) 


dt'. 


( 12 ) 


(13) 
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As the right hand side contains a convolution, taking the Laplace transform, 
£, of the equation with respect to time and rearranging gives, 


dKMi-c 


m 


1 


-C 


J o(t) 


e(t, o) j ~ 1 9{t, o) J c m)} \ o(t, o) j c m)} ■ 

This result can then be used to write, 

_ r -x (cm)} r i(t) 


(14) 




c { c{m} C l*M) 


_ C -1 f r \ MO 


Jo MM 0) 9(t',o)J 


cm)} \e(t, o) 


(15) 


where we have defined the memory kernel, 




Equation (12) allows us to write Eq. (11) as, 
dl(t ) 


dt 


= w(t)S{t)I{t)-6(t,0) 9 0 ^/^ dt'-'y{t)I(t)+0{t, 0) 


(16) 


d f I 0 (t) 


which, using Eq. (15), becomes, 


dt \6(t, 0) J 

(17) 


dl(t) 

dt 


= 0) K(t - t') 




0 (t',O) 9(t',0)J ” dt 

(18) 


M0 

9(t, 0) 


This equation is the generalised master equation that describes the time evo¬ 
lution of the number of infected individuals in an SIR model with arbitrary 
waiting time in the infectious compartment. 

Simple flux balance considerations give the master equations for the other 
two states. The equations for the susceptible and recovered populations are, 


dS(t) 

dt 

dR(t) 

dt 


= A(t) - u(t)S(t)I(t) - 7(t)S{t), 

/ m mo \ J4/ d (i 0 (t) 




lit 

0) 6(t', 0) J dt \0(t, 0) 


(19) 

-7 {t)R(t). 

( 20 ) 


Equations (18), (19), and (20) are the governing, or generalised master, equa¬ 


tions for an SIR model with a general recovery probability. 
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2.1 Relation to the Classic SIR Model 


The master equations for the SIR model with a general recovery probabil¬ 
ity reduce to the classic SIR model equations, with births and deaths, if 
the probability of not clearing an infection is exponentially distributed, i.e. 
</)(t) = exp(—/it). In this case the probability of an individual clearing an in¬ 
fection is not dependent on the amount of time that the person has already 
been infected. Subsituting the exponential distribution for <f>(t) into the kernel, 
Eq. (16), we obtain, 

K(t — t') = /i5{t — t'), (21) 

where S(t) is the Dirac delta function. Also noting that as <j>(t) = exp (—/it) we 
can write, 

± fy$_\ = _ hit) 

dt \9(t, 0 )) ^ 6(t, 0) 


( 22 ) 


We can now substitute the expression for the kernel, Eq. (21), into the gen¬ 


eralised master equations, Eqs. (18), (19), and (20), to yield the classic SIR 
equations, 


dS(t) 

dt 

dl(t) 

dt 

dR(t ) 
dt 


= A (t) - uj(t)S(t)I(t) - 7 (t) 5 '(t), 
= oj{t)S{t)I{t) - nI(t ) - 7 (t)I(t), 
= »I(t) 


(23) 

(24) 

(25) 


2.2 Relation to the Kermack and McKendrick Age-Structured Model 


The master equations for the SIR model with a general recovery probability 
are formally equivalent to a reduction of the general SIR model presented 
by Kermack and McKendrick [25]. The derivation of the Kermack and McK¬ 
endrick model from our stochastic process is presented in Appendix A. Here 
we show how the master equations, Eqs. (18 1 , (19), (20), can be obtained from 
a reduction of the Kermack and McKendrick SIR model equations given by, 


dS 

dt 

dt 

dt 

dR 

dt 

m 


A- 

di 

/>oo 

/ uii(a,t)daS(t) — "fS(t), 

Jo 

(26) 

da 

/»oo 

= -P(a)i(a,t) — 7 't(a, t), 

(27) 

l 

/»oo 

f3(a)i(a, t)da — 7 R(t), 

(28) 

L 

i{a , t)da. 

(29) 


In this model i(a,t) is the number of individuals who are infected at time t 
and who have been infected since time t — a. The equivalence can be seen by 
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making the identification, 

i(a,t) = &(t,t — a)q + (I,t — a). 


(30) 


Then Eq. (29) is equivalent to Eq. Q , provided that the separability assump¬ 
tion, Eq. ([ 6 j), holds. If we assume that i(a,t) —> 0 as a —> oo. Integrating Eq. 
(271 with respect to a then gives 


d.1 

dt 


fOO 

— i(0,t) = — / /3(a)<j>(a)9(t,t — a)q + (/,t — a)da — 7 /. (31) 

J 0 


Identifying fd{a)(f>{a) = ip (a) and i(0,t) = q + (I,t), we can then split the inte¬ 
gral to, 

dl f* 

— = i(0,t)— / ilj(a)0(t,t—a)q + (l,t—a)da— / ip(a)0(t,t—a)q + (I,t—a)da—'yl. 

dt Jo J t 

132) 

(18), hence with a change of variable of integration t' = t — a becoming. 

rt 


This allows for the use of the same Laplace transform method as in Eqs. (|14|)- 

able of integrat 

m w) 




0 ) 0 (^, 0 ) 


dt'+6{t, 0)^- ( 


dt \6(t, 0 ) 


Using the result that, 


««, o)4 ' m 


dt \0(t, 0 ) 


r° 

■I ij)(t — t')6(t,t')q + (I,t')dt'. 

J —OO 


- 7 1■ 

(33) 

(34) 


Hence we have recovered the generalised master equation given in Eq. (18) 
with time independent rates. 


2.3 Integral Equation Formulation 

The master equations for the fractional recovery SIR model can be formulated 
as a coupled set of integral equations. This enables comparisons with other 
related models and it enables the application of integral equation methods for 
analysis of equilibrium states. To formulate the system as integral equations 
we begin by noting that Eqs. 0 . @ can be substituted into Eq.Q to yield, 

q + {d,t) =u(t)S{t)I(t), (35) 

and then Eq. Q can be re-written to obtain the integral equation for the time 
evolution of the infected state, 

J(t)=I 0 (i)+ f $(t,t')u{t')S(t')I(t')dt'. (36) 

Jo 
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The integral equation for the time evolution of the susceptible state can 


be obtained by direct integration of Eq. (19), to yield, 

S(t) = 5(0) + f A(f') dt' - f u)(t')S(t')I(t') dt’ - f 7 (t')S(t') dt'. (37) 

Jo Jo Jo 


Note that the total population, 


N(t) = S(t) + I(t) + R(t), 


(38) 


and the master equations with a fractional order recovery were obtained with, 

d ^ = \ (t)-y(t)N(t). (39) 


We can combine Eqs. (38) and (39) to obtain the differential equation for the 


time evolution of the recovery state in the form, 


^jr = ^ + x (t) - 7 (t)S(t) - 7 - 7 (t)R{t), (40) 


and then integrate to find 


R(t) = R( 0) - lit) + 1(0) - 5(f) + 5(0) + f A(f') dt 1 

Jo 

[ 7(t , )5(t / ) dt 1 - f ^(t')I(t') dt'- f 7 (t')R(t') dt'. (41) 

Jo Jo Jo 

After substituting for I(t) and S(t) using Eqs.(36),(37) we have the integral 
equation for the time evolution of the recovered state, 

R(t) = R(0) + I 0 (0) - I 0 (t) - f <P(t,t')u(t?)S(t')I(t')dt' + f uj(t')S(t')I(t') dt' 

Jo Jo 

- f 7 (t')I(t') dt 1 - f 7 (t')R(t') dt 1 . (42) 

Jo Jo 


Equations (371, (36 1 , (421, provide a general set of coupled integral equa¬ 


tions for fractional recovery SIR models. The integral equation for the suscep¬ 


tible state, Eq. (47), can be shown to be equivalent to the integral equation ob¬ 


tained from S(t) = N(t) — I(t) — R(t) in the special case where A (t) = 7 (t)N(t), 
and thus N(t) = N( 0) is constant. 













10 


C. N. Angstmann et al. 


2.4 Reduction to the Hethcote and Tudor Endemic Disease Model 


If 7 (t),u(t) and A(f) are constant in time then the integral equations, (371, 
(36), ( |42| ) reduce to the integral equation model for endemic infection diseases 
that was introduced by Hethcote and Tudor J 2 T]. 

We first note that, with 7 constant, 


and, 




dY(t) d 


dt 


+ 7 Y{t) = e-^-{e^Y{t)) 


Substituting Eq.(43) into Eq.(36), with lj also constant, we have, 

I(t) = I 0 (t) + f (j)(t - dt'. 

Jo 


Using Eq. (44) we can re-write Eq. (19), with a;, 7 and A constant as, 

4 (e 7 *S(t)) = e 7t A - e 7 t wS(f)/(t), 

and then integrate with respect to time to obtain, 

S(t) = e _ 7 *S'(0) + f e-^-^Xdt' — f e- 7(^ - ^ ' ) a;5(^ , )/(^ , ) dt'. 
Jo Jo 


Using Eq. (44) we can re-write Eq. (40), with a;, 7 and A constant as, 
4 (e 7 t i?(t)) = Ae 7 * — 4 ( e 7 t 5 ( f)) — 4 ( e 7 *J(t)), 


(43) 

(44) 

(45) 

(46) 

(47) 

(48) 


and then integrate with respect to time to obtain, 

R(t) = R( 0)e" 7 * - I(t) + /(0)e" 7t - S(t) + S( 0)e" 7t + f Ae" 7 ^-*') dt'. (49) 

Jo 


We now substitute for I(t) and S(t) using Eqs.(45), (47) to obtain, 

R(t) = i?( 0 )e- 7 t +/ o ( 0 )e- 7 t -/ o (f)+ f ojS(t')I{t')e-^ t - t ''> (1 - </>(f-t')) dt'. 

Jo 

(50) 

Equations (451 and (50) recover the integral equations for the infected state 
and the recovery state in the endemic infection diseases model introduced in 
j2X|. The integral equation for the susceptible state in this case can be shown 
to be equivalent to the integral equation obtained from S(t) = N — I(t)— R(t) 
with N(t) = N( 0) = A/ 7 . 
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3 Fractional Recovery SIR Model 


When a person is persistently infected for a long period period of time, with 
little chance of spontaneous recovery, they are said to be chronically infected. 
This type of behaviour is not captured by the assumptions of the standard 
SIR model, i.e., exponentially distributed waiting times. We can incorporate 
chronic infections by having the probability of clearing the disease decrease 
with the amount of time that an individual has been infected. In a power-law 
tailed waiting time distribution the expected waiting time diverges. Using such 
a distribution in our SIR model will lead to individuals becoming “trapped” 
in the infectious compartment until they die. By utilising a power-law tailed 
Mittag-Leffler waiting time distribution our general SIR model will reduce to 
a set of differential equations with a fractional order time derivative on the 
recovery transition. 

The Mittag-Leffler probability density is defined by [23] . 

^(*) = (- (*) )> ( 51 ) 


with 0 < a < 1. Here E 0 ^(z) is the two parameter Mittag-Leffler function, 
defined by, 


E a,p{z) 


S r ( fca + ^)' 


(52) 


The Mittag-Leffler distribution limits to an exponential distribution in the 
case a = 1. For 0 < a < 1 the density has a power-law tail at long times ma. 


ip(t) ~ t~ x ~ a . 


(53) 


The corresponding survival function is given by, 


(t>[t) = E aA 



(54) 


The Laplace transform from t to s of the memory kernel, Eq. (16), for a 
Mittag-Leffler distribution is given by, 


C t [K(t)\s} = s'-V. 


(55) 


Here we use the notation £ t \Y(t)\s] to denote the Laplace transform of Y(t) 
from t to s and we use the notation £j 1 [F(s)|t] to denote the inverse Laplace 
transform of Y (s) from s to t. Without loss of generality we define, 


H = r“. (56) 

We also note that the Riemann-Louiville fractional derivative, defined by [3J 


Q V\~ a f{t) = 


m 


i 


r(a) dt J 0 (■t — t') 1- 


dt' 


(57) 
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has the Laplace transform 


A[ o^t 1 a /W|s] = s 1 “A[/(i)|s] 


(58) 


so that the fractional derivative can also be written as the following convolu¬ 
tion, 


(Pi 


7W = / 

Jo 


C 


-i r s i-«u/ 


t 7 ] f(t — t')dt'. 


(59) 


It follows from Eq. (55) and Eq. (59) that if the kernel in Eqs. (18), (19), (20) 
is obtained from the Mittag-Leffler waiting time density then we obtain the 
generalised master equations with fractional recovery, 


dS{t) 

dt 

dljt) 
dt 

dR(t) 
dt 


= A (t) - u(t)S(t)I(t) - 7 (t)S(t), 

= w(t)S(t)I(t) - 7 (t)I(t) - 9(t, 0) ( fiD. 


I(t) hit) 


= e(t, o)( / xP t 1 -“ 


i(t) hit) 


9(t, 0) 9(t, 0) J dt \9(t, 0) 


9(t, 0) 9{t, 0) 

hit) 


(60) 

/ ^o(f) 
dt \9(t, 0) 
(61) 


-7(f)#(f)- 


(62) 


In the following we will refer to the above set of equations as the frSIR model. 
Letting a = 1 recovers the standard SIR model. 


3.1 Equilibrium States 


The frSIR model is a non-autonomous dynamical system which can be simpli¬ 
fied, by taking the birth, infectivity and death rate to be constant parameters, 
A {t) = A, ui(t) = ui and 7 (t) = 7 respectively, giving, 


dS ^ = A- wS(t)I{t) -'ySit), (63) 

^ = cosmt) 7 m - e_7< [doV]~ a (e 7t (7(i) - I 0 (i))) - | (e 7t / 0 (f))) , 

(64) 

^ (e 7 ‘ (J(i) - Jo(t))) - j t (e^Io(t))^ - 7 R(t). (65) 

We take the equilibrium state to be iS*,I*,R*), such that, 

lim Sit) = S*, lim I it) = /*, lim R(t) = R*. (66) 


Taking the limit as t —» 00 of Eqs. (63), (64), and (65), and noting that, 
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we have 


0 = \-uS*I* - 7 S*, 

0 = uS*I* - lim e" 7 Vo^ t 1_ “ (e 7 ‘ ( lit) - I 0 (t ))) 
0 = lim e - 7 Vo2? f 1_a (e 7 ‘ (/(t) - J 0 (t))) - lR* ■ 

£—>■ OO 


( 68 ) 

(69) 

(70) 


In order to calculate the unevaluated limits in Eqs. (691, and (70) we 
consider a Laplace transform of the terms, 


£{e- 7 t o^-“ (e 7t J(f))} = (s + 7 ) 1 - 0 (•/(*)) • 
In which we have defined, 


(71) 

(72) 


J(t) = I(t) - I 0 (t). 

We can express Eq. using a Taylor series expansion, 

J(s)(s + 7 y~ a = J(s) (7 1_a + (1 - a) 7 - a s + 0(s 2 )) . (73) 

As the Laplace transform is a linear operator we can take the inverse termwise, 
producing, 

e" 7t o V\~ a (e 7t (/(t) - I 0 {t))) = C~ l {J(s) ( 7 l ~ a + (1 - a) 7 ~ a s + 0(s 2 ))}, 

(74) 

= 7 1 "“J(t) + (1 - a) 7 ~ adj P- + C - 1 ( 0 (s 2 )) . 


dt 


The limit of J(t) is, 


lim J(t) = I*. 

£—>■00 


(75) 


(76) 


It is then clear that, in the long time limit, the inverse Laplace transform of 
the higher order terms of the Taylor expansion will become zero, i.e. 


lim £ 2 ! = 0 , 
*-►0 dt 

lim ( 0 (s 2 )) = 0 . 
t-K> v ' 


Thus we can compute the desired limit, 


lim e- 7 ‘ o V\~ a (e 7t (I(t) I 0 (t ))) = ^1* 

£—>■00 v 7 


Substituting Eq. (79) into Eqs. (|63j), (64), and (|65j) yields, 

0 = A - ojS*I* - 7 S*, 

0 = ojS*I* — /z 7 1 _a J* — 7 /*, 

0 = iij 1 -*!* - 7 R*. 


(77) 

(78) 

(79) 


(80) 

(81) 

(82) 
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Solving Eqs. (80), (81), and (82) reveals two equilibrium states, the disease 
free state, 


S* = ~, I* = 0 , R* = 0 , 


(83) 


and the endemic state, 


y = W*r"+T , j, = *-2. -gT! (84) 

w /xq 1 a + 7 w /X7 + 7 1+a w 

As the population in each compartment can not be negative the endemic 
equilibrium can only exist if, 


Xoj > /Z 7 2 “ + 7 2 . 


(85) 


When a = 1, the equilibrium states are equivalent to the fixed points of the 
classic SIR model with constant parameters. For 0 < a < 1 the equilibrium 
states are not fixed points. This invalidates the use of a standard linear stability 
analysis around the equilibrium states. However some progress can be made 
by considering the integral equation formulation of the frSIR model, Eqs. (451, 
(471, and (50). After a translation of the equilibrium states to the origin, the 
asymptotic behaviour of the resulting system of nonlinear Volterra integral 
equations is equivalently given by the asymptotic behaviour of its linearization 
as a system of linear Volterra integral equations JIT| ■ This remarkable result 
was used by Hethcote and Tudor [31. to infer local stability properties of the 
equilibrium states of the integral equations Eqs. (45), (47), and (50). If the 
results of Hethcote and Tudor [21] are applied to the special case of the frSIR 
model, where </>(t) is defined in Eq.(54), then the disease free state, Eq.(83), is 
locally stable if Aw < n r y 2 ~ a + 7 2 and the endemic equilibrium state, Eq.( 84), 
is locally stable when it exists. 


4 SIR as a Discrete Time Random Walk 

There are numerous numerical methods that have been developed for solving 
fractional order differential equations [33113511151 and many of these methods 
could be adapted to the system under consideration here. However recently 
we showed that in the case where the fractional order derivatives have been 
derived from continuous time random walks it is useful to reformulate the 
problem using discrete time random walks (DTRWs) and then use this for¬ 
mulation as the basis of a numerical method [3]. The advantage of basing the 
numerical method on a discrete time stochastic process is that the derived 
explicit numerical method is easy to implement and is also inherently stable. 
























A fractional order recovery SIR model from a stochastic process 


15 


4.1 Discrete Time Random Walk 

We consider a discrete time random walk where the walking particle is an 
individual who will transition though the S, I, and R, compartments. In order 
to obtain a useful numerical scheme from this we need the discrete time process 
to limit to the continuum process as the time step, At, goes to zero. This 
necessitates a slight modification to the transitions that were considered in 
the continuum case. Individuals are born into the susceptible compartment 
with a birth rate A (t) so that the number of individuals being born on the 
n th time step, between time t and t + At, is equal to A(n) = A (t)At. The 
probability that an individual will die and be removed from consideration on 
the n th time step is 7 (n). Susceptible individuals who come in contact with 
an infected individual may become infected. The probability of an infected 
individual coming in contact and infecting a susceptible individual in the n th 
time step is u>(n). 

The recovery of infected individuals is assumed to be dependent on the 
number of time steps since they entered the infectious compartment. The in¬ 
dividual is assumed to wait in the infectious compartment with a probability 
of “jumping” that is dependent on the time step. When a transition event oc¬ 
curs there are two possible things that can happen to the individual, they will 
either move to the recovered compartment or re-enter the infectious compart¬ 
ment with the transition probability, dependent on the time step, being reset. 
If the individual was infected before the first time step they will always tran¬ 
sition to the recovered compartment. This self jump modification is required 
to ensure appropriate scaling as the size of the time step tends to zero. The 
probability of transitioning to the R compartment, given a jump occurred, is 
denoted by r. The probability of jumping on the 71 th step, conditional on not 
having jumped in the first (n — 1) steps, is denoted by /z(n). 

From this it follows that the probability flux entering the I compartment 
on the n th time step can be written as 

n—1 

Qj(?z)= u>(n)S(n — l)(j)(n — 1 — k)0(n — l,k)Qf(k) 

k =-°° ( 86 ) 

+ (!-»•)£ /r(n — k)4>(n — 1 — k)9(n , k)Qf ( k ). 

fc=1 


In the above equation, <p(n — k) is the probability of not jumping for (n — 
k) steps and 0(n, k) is the probability that an individual who entered the I 
compartment on the k th step has survived the death process up to the n th 
step. If we assume that we have an initial distribution of infectious individuals 
at the 0 th time step, then we can infer that the flux at earlier times is given 

by, 


Q+(n) 


i(—n, 0 ) 
9(0, —n)<j>(—n) 


Vn < 0. 


(87) 
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Here i(n, 0) is the number of individuals at time step 0 who have been infected 
since time step n , with n < 0. This allows us to write the flux for n > 0 as 


Qf (n) =Lu{n)S{n — 1) </>(n — 1 — k)9(n — 1, k)Qf (k ) 

fc=1 

n— 1 

+ (1 — r) /r(n — k)<p(n — 1 — k)9{n , (fc) 


( 88 ) 


fc=i 


+ u)(n)S(n — l)0(n — 1, 0) ^ (j)(n - 1 - 


k=—c 


We can easily see that 


i—k 


(j>{n — k) = J|(l - mC0)> 


(89) 


j=0 


and 


0{n, k) = n (! ~ TO'))- (9°) 

j=k 

The number of individuals in the / compartment on the n th time step is 


where, 


7 ( n ) =£ (j)(n — k)6{n , ( k ) + /o(n). 

fc=i 

I 0 (n) = 6»(?r,0) ^ <j>(n~ k) l ^ k ^ . 


k ——oo 


(91) 


(92) 


Subtracting J(n — 1) from each side of Eq. (911 gives, 

n— 1 

I(n) — J(n — 1) = Q/(n) + ^^(0(n — k)9(n, k) — </>(n — 1 — k)9(n — 1, k))Q~\ ( k) 


k =1 


+ / 0 (n) - / 0 (n - 1), 


(93) 


and substituting Eq. (88) into the right hand side of Eq. (93) and using Eq. (91) 
gives, 

I(n) — I(n — 1 ) =oj(n)S(n — 1 )I{n — 1 )) + Io(n) — Iq{ji — 1 ) 

n— 1 


+ (!-»•)£ /r(?z — k)<j>(n — 1 — k)9{n , k)Qf (k) 


k =1 


+ — k)9{n , fc) — <j)(n — 1 — k)9(n — 1, k))Qf ( k ) 

( 94 ) 


fc=i 
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Noting that, 

(4>(n—k)9(n, k)—<j>(n—l—k)6(n—l, k )) = — r y(n)(f>(n—l—k)9(n—l, k)—n(n— k)<j>(n— 1— k)9(n, k). 

(95) 

we get, 

I(n) — I(n — 1) =uj(n)S(n — 1 )/(n — 1) + Io(n) — Io(n — 1) — 7 (n)(/(n — 1) — /o(n — 1)) 

n— 1 


— r fj,(n — k)(j>(n — 1 — k)9(n, k)Qf (fc). 


fe=i 


(96) 


To simplify further we mirror the derivation in the CTRW approach, using 
the semi-group property of the death survival function, and using discrete Z 
transform methods to replace the memory kernels. The Z transform form n 
to z of Y (n) is defined by 


Z[Y(n)\z] = J2 Y (n)z- n . 


n —0 


First we use the result 


6(n, 0) = 9(n, k)9(k , 0), 


(97) 


(98) 


in Eq. (96) to write 


J(n) — J(n — 1) =u>(n)S(n — 1 )I{n — 1) + Io(n) — /o(n — 1) — 7 (n)(/(n — 1) — Io(n — 1 )) 


0 + (k) 

- r9(n, 0) Y »(n - k)(j){n - 1 - k) ■ 


k=1 


and we use the same result in Eq. (91) to write 
I(n) - 7 0 (n) 


9{n 1 0) 


X! i *9/ W 

= E^- fc )^y- 


(99) 


( 100 ) 


k=1 


Finally, we need to express the sum n{n—k)(j){n— 1 —fc) , in terms of 

I(n) so we take the Z-transform of Eq. (100), and use the convolution property 
of Z transforms to give, 


Ijn) - I 0 (n) 
9(n, 0) 1 


= Z [<j){n)\ z\Z 


We can now write 

n— 1 


Y n - - 1 - = Y K ( n ~ 1 - 


r Q/(n) 

9(n, 0) 

I(k) - I 0 (k) 


( 101 ) 


k=0 


0(fc,O) 


fc =0 


9(k, 0) 


( 102 ) 
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where, 


c (n) = Z 


Z [/j,(n)cp(n - l)|z] 


(103) 


Z[<t>{n)\z] 

The discrete time generalised master equation for an arbitrary waiting time 
distribution can then finally be found by substituting Eq. (1021 into Eq. (99) 
to give, 

I{n) — I{n— 1) =ui(n)S(n — 1 )I(n — 1) + Io(n) — /o(n — 1) — 7 (n)(/(n — 1) — 


m nl - Ip(k) 

-rf(n,0 )£*(„-k) - 


k=1 


(104) 


The master equations for the other compartments can be found by again con¬ 
sidering a flux balance so that, 


S(n) — S(n 
R(n) — R(n 


1) = A(n) — uj(n)S(n — 1 )I(n — 1) — "f(n)S(n — 1), 


k =1 


+ I 0 (n) - I 0 {n - 1) + 7(n)/ 0 (n - 1) - 7 (ri)R(n - 


(105) 

(106) 
!)• 


4.2 Discrete Time Fractional Recovery SIR Model 

The continuous time frSIR model was obtained by considering Mittag-Leffler 
waiting time densities in the CTRW formulation. In the DTRW formulation 
a discrete time fractional recovery SIR model is obtained by considering a 
Sibuya(a) waiting time distribution [39ll4], In this case the probability of jump¬ 
ing on the n th time step, conditional on not having jumped on the previous 
n — 1 time steps, is given by 


M«) = 


0, n = 0, 

and the jump survival probability is given by 

T(1 — a + n) 


’Pin) = 


r(n + i)r(i - a)' 

It is a simple matter to obtain the Z transforms from n to z 

a -1 


Z [<p(n)\z] = 


z - 1 


(107) 


(108) 


(109) 


I 0 {n - 1)) 
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and 


Now using Eq. p3| we have 


^2 H(n)<l>(n - l)z n 

(110) 

n—0 


+(n)0(n - l)z~ n 

cm) 

(z- 1\“ 


*-( . ) ' 

(112) 


(113) 


r—1 


and then after taking the inverse Z transform we obtain the fractional memory 
kernel 

, x r , r(n - 1 + a) 

Kl - n)=i '-~ + na-mn + iy (114) 

This memory kernel can be calculated recursively by noting that for n > 3, 

a — 1 


k(ti) =11 + 


n{n — 1). 


(115) 


The first two values are simply k(2) = f(a — l)and k(1) = a. 

The number of initially infected individuals left at the n th time step, Io(n), 
is found from Eq. (92) for the case 0 < a < 1, 

, . . +(1 — k)r(l — k + n — a) , . . 

Jo(n) = »(n,0) E r(l _ k + n)r{l _ t - < 116 ) 

In the case where the initially infected individuals were all infected at time 
zero this will simplify. Taking 0) = ioSk,o where io is a constant and 5k,o 
is a Kronecker Delta function, we find, 


r / \ m r(l + n-a) . 

I 0 {n) = 0(n, 0) -rz 0 . 

1(1 + n)I (1 — a) 

Again this may be also expressed recursively, 

Io(n) = (l - Io(n - 1), 


(117) 


(118) 


with the initial condition, /o(0) = Iq. 

When a = 1, <j>(n) = 0, for all n > 0, as such care has to be taken with the 
definition of Io(n). In this case the only physically permitted initial condition 
is that the initially infected individuals were all infected at time zero. Hence, 
when a = 1, Io{n) = 0 for n > 0, and lo(0) = io- 

The discrete time fractional recovery SIR model is obtained by using the 
memory kernel, Eq. (114), in Eqs.(104), (105), and (106). 
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In Appendix B we show that the discrete time fractional recovery SIR 
model equations limit to the frSIR model equations by identifying t = nAt 
and taking the limit At —>■ 0 and r —>■ 0, with 


lim —— = u. 

zlt,r->-0 At a 


(119) 


This justifies our use of the explicit discrete time model equations as a numer¬ 
ical method for solving the continuous time frSIR model. 


4.3 Equivalence Between The Discrete and Continuous Parameters 


As the discrete frSIR governing equations, Eqs. (105), (1041, and (106) with 


Eq. (114), limit to the continuous time frSIR model equations we can approx¬ 


imate the solution of the continuous time equations with the solution to the 
discrete time equations. Given a set of parameters for the continuous time 
equations we need to identify corresponding values of the parameters in the 
discrete time equations. The continuous time frSIR model is parameterised 
by three functions, (A(t), u>(t), 7(f)), and two constants ( a , /x). The discrete 
time frSIR model is parameterised by three functions (A(n), ui(n), 7(71)), and 
three constants (a, r, At). Assuming a given time step At, the correspondence 
between the continuous time t and the discrete time n is given by t = nAt. 
The expected number of births in a time step is related to the continuous time 
birth rate by, 

A(n) = AtX(n5t). (120) 

The probability of an individual becoming infected in a time step is related to 
the continuous time infection rate by, 


/ Ann 

;(n) = 1 — exp — / 

y J nSt 


( n-\-l)At 


;(0 dt' 


( 121 ) 


Note that in the limit At —> 0 Eq. (121) is consistent with Eq. (156). We can 


treat the death probability in a similar fashion, 


/ j.(n+l)At \ 

7 (n) = 1 - exp - / 7(0 dt' . 

\ JnSt ) 


( 122 ) 


This definition of the death probability ensures that the discrete survival func¬ 
tion, 6{n , m) is equal to the continuous time survival function evaluated at nSt 
mSt, i.e. 


( r-nAt \ 

— / 7(0 dt' ) . 

J mAt J 


(123) 


The discrete anomalous exponent, a, is unchanged from the continuous case. 
The remaining parameter in the discrete model, r, is obtained from 


r = /j,At a 


(124) 
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which is consistent with Eq.(119l. Note that r is a probability and as such 


is bound in the interval [0,1]. Thus given the parameters and a from the 
continuous time model, Eq. pi) constrains At as follows, 


At< ( — 


(125) 


Given that Eqs. ( |105 ), ( 104| ), and (106) are explicit difference equations if 
we know S(k), I{k), and R(k), for 1 < k < n, then we can calculate S(n), 
I(n ), and R(n). In this manner we have a simple numerical method that 
approximates the solution of the continuous time frSIR model by simply noting 
that for t = nAt, S(t) ss S(n), I(t) ~ /(n), and R(t) ss R(n). 


5 Example 


We consider the frSIR model, Eqs. (60), (61), and (62), with the following 
parameters, A (t) = 0.1, uj(t) = 0.02, y(t) = 0.001, n = 1, and a range of 
a £ (0,1]. We also take a delta function initial condition at t = 0 for the 
infected population, i(t, 0) = 0.5<5(f). The initial recovered population is taken 
to be zero, -R(O) = 0, and the initial population in the susceptible compartment 
is taken so that the total population is at the equilibrium level of 100, i.e. 
S( 0) = 99.5. With these parameter values, the system has two possible steady 


states, the disease free steady state given by Eq. (831 and the endemic steady 


state given by Eq. (84). The endemic steady state for susceptibles and infecteds 
are plotted as a function of a in Fig.Q. 

To find the numerical approximation for this situation we solve the dis¬ 
crete equations, Eqs.( 105| ) , (104), and ( |106 ). For a given At, the discrete pa¬ 
rameters are taken to be A(n) = O.lAt, oj(n) = 1 — exp(—0.02Z\f), y(n) = 
1 — exp(—0.001Z\f), and r = At a . The initial conditions are implemented by 
taking 5(0) = 99.5, I?(0) = 0, and *(— k, 0) = 0.5<5o,fc- 

In Fig. Q we show plots of 5(f) and /(t) versus time for a = 0.3,0.5, 0.7,0.9 
with initial conditions 5(0) = 99.5,/(0) = 0.5,i?(0) = 0. The number of sus¬ 
ceptibles falls sharply over short times before rising slowly towards a steady 
state at later times. The number of infecteds rises to a maximum in the short 
time regime before declining slowly back towards a long time steady state. It 
can be seen that the peak level of infection and the long time levels of infection 
are both increased with decreasing a. 

For values of a very close to one the numerical solutions show an oscillatory 
approach towards the steady state. This can be seen in Fig. © for a = 0.99, 
0.999, and 1. The oscillations are suppressed as a is decreased. 


6 Summary 

We have derived a fractional order recovery SIR model that differs from the 
classic SIR model by considering the case where the recovery from a disease 
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a 

Fig. 1 The endemic steady state in the frSIR model plotted as a function of a for sus- 
ceptibles (solid line) and infecteds (dashed line). The Parameters were A = 0.1, to = 0.02, 
7 = 0 . 001 . 






Fig. 2 Plots of S(t) and /(£) versus time in the fractional recovery SIR model with a = 
0.3, 0.5, 0.7, 0.9. The arrow indicates the direction of increasing a. The other parameters 
are A = 0.1, co = 0.02, 7 = 0.001. The model was solved using the DTRW method with 
At = 0.05. 
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Fig. 3 Plot of susceptible versus time, after the initial relaxation time, for a = 1 (solid 
line), a = 0.999 (dotted line) and a = 0.99 (dashed line). The other parameters are A = 0.1, 
u = 0.02, 7 = 0.001, and At = 0.05. 


does not follow an exponential distribution, but follows a distribution with 
a power law tail. The fractional order model is biologically motivated by the 
observation that in some disease processes, the longer a person is infectious 
the more likely they will remain infectious. The fractional order model per¬ 
mits both a disease free equilibrium state and an endemic equilibrium state. 
We have related the fractional order model to the generalised SIR models in¬ 
troduced by Kermack and McKendrick. The fractional order model that we 
have derived is different to ad-hoc fractional order epidemic models with Ca- 
puto fractional derivatives in time replacing standard derivatives in time. Our 
work justifies an appropriate inclusion of fractional derivatives in SIR model 
equations. 


We have also derived a discrete time fractional order recovery model that 
limits to the continuous time fractional order model as the time step goes to 
zero. The discrete time model was itself derived from a stochastic process and 
we have used this model to provide a stable numerical solver for the continuous 
time model. Our numerical solutions based on this discrete model show that 
the prevalence of infection and the peak levels of infection are both elevated 
by the fractional order derivative as it is varied further away from the integer 
order derivative in the classic SIR model. 
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Appendix A 

Derivation of the Kermack-McKendrick Model for Infecteds in a 
Age-Structured System 


Let i(a, t) denote the number of individuals who are infected at time t and 
who have been infected for time a. The total number of infected individuals 
at time t is given by 


; ft 

i(a,t)da= / ^(i, t') dt' 

Jo 


(126) 


where 


0(t,l/) = 4>(t-t')e{t,t') (127) 

is the product of the probability of surviving death, and of ‘surviving’ 

the transition out of the next stage, <p(t — t'). If we differentiate Eq.( 126) with 
respect to t using Leibniz rule we have 


l °i da+iM= l 


d ^ t: 0, t') dt' + $>(t, t)i(0, t) (128) 


dt 


We also have 


But 


f* di 

— da = i(t , t) — i(0, t). 


>o 


<P(t,t) = 1 


so that if we add Eq.(128) and Eq.(129) then we obtain 


/o 


* / di di , 

m + ei da = 


0 


dt 


i( 0, t') dt' 


(129) 

(130) 

(131) 


We can expand the partial derivative of $(t,t r ) using the product rule in 
Eq. (127) with 

d0(t, t') 


and 


to obtain 


dt 

d(j>{t — t') 
dt 


= -7 (t)0(t,t r ) 
= -ip{t-t'), 


(132) 

(133) 


— + —da = — f ip(t—t')6(t,t')i(0,t')dt'—f 'y(t)6{t,t')<f>{t—t')i{Q,t')dt'. 
at oaj J o J o 

(134) 

Using Eq.( 1271 and then Eq.( 126) in the second term on the RHS of Eq.( |134 ) 
we now have 

J ^ da = — J i/)(t—t')6(t,t')i(0,t')dt'—'y(t) J i(a,t)da. (135) 
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We now define 


P(t-t') = 


y >(* - ?) 

- 1 ') 


(136) 


and use Eq.( 1271 in the first term on the RHS of Eq.( 1351 to obtain 


— +—} da = — I j i(a,t)da. (137) 


In the first integral on the RHS of Eq.(137) it is convenient to make a change 
of variables t' = t — a with dt' = —da, then we have 


di di 

— + — ) da = — j (3(a)d>(t,t — a)i(0,t — a) da — 7(f) / i(a,t)da. 


Finally we note that 


i(a, t) = <I>(t , t — a)z(0, t — a) 


so that Eq. (138) can be written as 


(138) 

(139) 


— + —^ da = — f /3(a)i(a,t) da — 7(f) [ i(a,t)da. (140) 
at da In Jn 


This is consistent with Eq. (27) for the change in infectives in the age-structured 
Kermack McKendrick model. 


Appendix B 

Limits To Continuous Time 

The discrete time fractional recovery SIR model can be shown to limit to 
the frSIR model by identifying t = nAt and taking the limit At 0 with 
r/At a finite. The continuous time equations can be obtained from the discrete 
time equations using Z star transform methods. The Z star transform of Y(n) 
is given by 

OO 

Z*[Y(n)\s, At] = Y(n)e~ nsAt (141) 

n—0 


Y{n)e~ nsAt At 

n—0 

00 

(142) 

Y(nAt)e~ nsAt At 

n—0 

(143) 


It follows that 
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where we have introduced Y(t) as a function defined over a continuous variable 
t. We can now take the inverse Laplace transform from s to t 


C 


-l 


AtZ*[Y{n)\s,At\ 


= Y ( nAt)S(t — nAt)At 


(144) 


n—0 


where S(t) is the Dirac delta function. Here, and in the following, we use the 


notation C 


-l 


Y(s) 


to denote the inverse Laplace transform from s to t and 


we use the notation C* 


Y(t) 


to denote the Laplace transform from t to s. 


It is useful to define the function 

OO 

Y(t\At) = Y ( nAt,)S(t — nAt)At. 

In a similar fashion we have 

C- 1 AtZ*[Y{n-l)\s,At] 


(145) 


n —0 


= y; Y ( nAt)S(t — (n + 1 )At)At (146) 

n—0 

oo 

= y; Y ((n — 1 )At)5(t — nAt)At (147) 


n=0 

= Y(t- At\At). 


Note that, with t’ = nAt , in Eq.( 145), we have 


lim Y(t\At) — lim V'' Y(nAt)S(t — nAt)At 

At ->-0 


At —>-0 


n—0 


Y(t')8(t - t') dt' 


= Y(t). 


This formally identifies 


Y (t) = lim C 
At->0 


AtZ*[Y(n)\s,At\ 


(148) 

(149) 

(150) 

(151) 

(152) 


provided that the limit exists. 

We further note the product rule 

OO 

^lim X(nAt)Y(nAt)8(t — nAt)At 


n—0 


= | lim X(nAt)5(t — nAt)At j ( lim Y(nAt)5(t — nAt)At\§ 3) 

\ n—0 ) \ ~~n=0 ) 
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which equates to X(t)Y(t) in each case, with t' = nAt , provided that both 
X(t) and Y(t) exist. 

We n ow take the inverse Laplace transform of the Z star transform of 
Eq.(104) and multiply by to write 


E 

n—0 


I(nAt) — I((n — l)Z\f) 


At 


S(t — nAt) At 


=E 


n—0 


Co(nAt) 

At 


S((n — 1 )At)I((n — l)At)5(t — nAt)At 


- ~^ t E o( nAt i °) (E R H n - k ) At ) 


I(kAt) - I 0 (kSt) 


n =0 


\k =0 


-E 1 ^ (/((n-l)Z\i)-/o((n- 


0(&Z\i,O) 

<5(f — nAt) At 


S(t — nAt) At 


n —0 
oo 


Y- I 0 (nAt) - / 0 ((n - 1 )At) 

+ 2^-- 6 ^ ~ 


n—0 


(154) 


We now take the continuous time limit of Eq.( 154) using t' = nAt and the 
product rule in Eq. (153), to obtain 

/»oo 

Jo 

— ( [ 9(t r , 0)5(t — t') dt , '\ ( lim EE I E H((n — k)At)~ 

vo ) *(***. °) 

- f (j{t') - io (t')) £(* -0 


/(fcZif) - Jo(fcZit) 


r°° / lim WEM 6(t-t’)dt’ 


/\t —>-0 


At 


where we have defined continuous time rate parameters 

Co (nAt) 


oo 


( t') = lim 


zit—>-o At 


and 


„ . 7 (nAt) 

7 (t') = lim 


(155) 

(156) 

(157) 
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Equation (155) simplifies further to 


dl(t) 

dt 


dlpjt) 

dt 


= u(t)S(t)i(t) - 7 (t) (i(t) - / 0 (t)) 

- i(t, 0) ( lim A £ (’g S((„ - S(t - nAt)At 


n —0 \k =0 


9{kAt, 0 ) J 


(158) 


The further reduction of this equation depends on the specific form of the 
memory kernel n{n). In the case of a jump at each time step the memory 
kernel is, 


K (n) = <5„, i. 


(159) 


In this case we can perform the sum over k explicitly in Eq. (1581 to arrive at 


df(t) 

dt 


= - 7 (t) ( fa ) - m) + ^ 

- 0) ( lim 1- V 

\ At ^ oAt ^p 0(( n — l)At, 0 ) 


^^-5(t-nAt)At) 


(\60) 




In order for the continuous time limit of the above equation to exist we define 


M = 


lim ——. 
zit->o At 


(161) 


Note that r is a free parameter in the range [0,1] and hence /i is only well 
defined in this limit if we take r to be a function of At. With this definition of 
/r we can now perform the limit (5i —>■ 0 to obtain the continuous time equation 


^ = u,(t)S(t)I(t) - M (f(t) - I 0 (t)) - 7 (t)(I(t) ~ Io(t)) + dI ^- (162) 

This further simplifies to 


^ = w(t)S(t)I(t) - nl{t) - 7 (t)i(t) 


(163) 


Equation (1631 recovers the corresponding equation in the classic SIR model. 


We now consider the continuous time limit of Eq. (158) with the Sibuya 


memory kernel, given by Eq.( 103). First we simplify the double sum in Eq.( 158) 
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using Laplace transforms and Z star transforms as follows: 

oo / n— 1 

At-* o At 


lim 1- V (V «((„ - <5(f - 

' 0(fc4t,O) y v ; 


n—0 \k—0 


= £ 


-l 


£t 


Km ^ V f V «((n - k)At) I{kM) h[kAt) ) 5{t nAt)At 
u ; 9{kAt,0) ) 


= £ 


= c: 


-i 


n —0 \k —0 
oo / n—1 


j“„ sE(E *«» - e- 


n—0 \k—0 
n— 1 




fc=0 


9(k, 0) 


9(kAt 1 0) y 
t 


snAt At 


-- £ -' J!”„ 5 z- W">' ! ’ |», a t ]a t 


The last line in the above follows from the convolution theorem for Z star 
transforms. 

To proceed further we use the Z transform of the Sibuya memory kernel in 


Eq.( 113) to write 

Z*[n(n)\s, At] = [((1 - e -^i-a _ (1 _ e ~sAt^ 


(sAt) 1 


j(sAt). 


(165) 

(166) 


The result in Eq.(164) can now be written as 


hm V ( V S((n - k)At) — T o( kAt ) ) 5{t _ nAt ^ At 
At^Q At A—t l .. U ; 9{kAt, 0) j 


n—0 \k—0 


-1 


= £. 
= ^s 1 

= i-i£s l 


lim f,i- y 7(n ^ } ~ 

At^oAt* ^ 9(nAt, 0) 

r-a r°° Hi) ~ Io(t) ^ st 


s C t 


9{t, 0) 

m - i 0 (t) 


-e dt 


9{t, 0) 


where 


u = lim ——. 
r At-vo Z\t“ 


(167) 


(168) 


Finally we substitute the result of Eq.(167) into Eq.( 158), and use the 
known result 1551 


(164) 


D l t ~ a Y(t) 


= 5 


1—a 


c* 


Y(t) 


(169) 
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to invert the Laplace transform and obtain 




m-m 
, S(t,0) , 


- 7 (t) (i(i)-I 0 (t))- 


dl 0 (t) 

dt 


(170) 

Equation (1701 recovers the continuous time frSIR model equation. 

Note that in order for the continuous time limit of the frSIR model equation 
to exist we defined 


u = lim —— 
^ At-vo At a 


(171) 


which requires r £ [0,1] to be a function of At. This is important for numerical 
simulations based on this DTRW method where we take r = )iAt a and then 
the requirement that r £ [0,1] places restrictions on At for given a and /i. 
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